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ABSTRACT 


An  efficient  machine  implementation  of  the  1972  COSPAR  International 
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Reference  Atmosphere  (CIRA)  is  developed  for  use  in  evaluating  atmospheric- 
drag  effects  upon  earth -orbiting  satellites.  Restriction  of  this  procedure 
to  the  satellite-orbit  region  of  significant  drag  and  the  use  of  an  alternate 
temperature  profile  representation  have  enabled  a  factor  of  30  improvement 
in  machine  execution  time  compared  with  that  required  by  the  code  supplied 
by  CIRA  publication  while  incurring  only  a  few  percent  difference  in  the 
density  results. 
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I.  INTRODUCTION 

This  note  documents  a  computer  code  entitled  CAD  for  the  Calculation 
of  Atmospheric  Density  at  satellite  altitudes.  This  code  has  been  developed 
specifically  for  use  by  the  ANODE  (Analytic  Orbit  Determination)  Program 
(Ref.  1) ,  which  is  currently  used  at  Millstone  Hill  as  a  rapid  orbit  esti¬ 
mator,  to  enable  estimation  of  atmospheric  drag  effects  upon  satellite 
orbits,  though  the  uses  of  this  code  could  be  more  general. 

CAD  is  intended  to  provide  an  efficient  machine  implementation  of  the 
CIRA  (Cospar  International  Reference  Atmosphere)  1972  (Ref.  2)  density 
model.  The  CIRA  publication  provides  a  computer  code  entitled  ADEN  for 
this  computation,  but  this  code  was  "intended  primarily  as  a  reference, 
rather  than  as  a  working  routine ,  ...  in  order  to  make  the  subroutine  as 
simple  and  as  closely  related  to  the  abstract  description  as  possible .  ... 

The  subroutine,  therefore,  runs  relatively  slowly  -  probably  too  slowly  for 
extensive  use  . . .  one  call  . . .  requires  the  equivalent  of  about  34000 
floating  point  multiplications"  (Ref.  2).  Most  of  this  calculation  time  is 
emended  in  the  numerical  integration  of  the  diffusion  equation  from  a 
lower  boundary  to  the  satellite  altitude.  CAD  replaces  this  numerical  inte¬ 
gration  with  an  analytic  approximation,  and  it  is  shown  that  the  difference 
between  the  two  results  is  probably  insubstantial  in  comparison  with  the 
uncertainty  in  such  atmospheric  models. 

Reference  1  defines  the  "significant  drag  regime"  of  the  atmosphere 
as  being  below  500  km  altitude.  The  CIRA  considers  hydrogen  to  contribute 
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negligibly  to  the  total  atmospheric  density  below  this  altitude,  and  thus 
the  portion  of  the  CIRA  dealing  with  hydrogen  may  be  omitted  for  the  present 
purposes.  Similarly,  if  we  adopt  125  km  as  the  lower  boundary  of  interest 
for  satellite  orbits,  a  not  unreasonable  choice,  then  we  may  avoid  the 
relatively  complex  portion  of  the  CIRA  applicable  to  altitudes  below  this 
level.  In  addition,  the  constituent  argon,  though  considered  in  the  CIRA, 
never  contributes  appreciably  to  total  density,  and  this  portion  of  the  CIRA 
may  also  be  omitted. 

Included  in  the  appendix  is  the  FORTRAN  code  for  procedure  CAD.  This 
code  is  a  straightforward  modification  of  that  given  in  the  CIRA  publication. 
Displacement  of  the  lower  altitude  boundary  from  90  to  125  km,  the  neglect 
of  hydrogen  and  argon,  and  the  substitution  of  an  analytical  calculation  for 
the  numerical  integration  of  the  diffusion  equation  are  the  only  changes 
made  to  the  published  CIRA  code. 

Section  II  of  this  report  discusses  the  problem  of  the  efficient  cal¬ 
culation  of  the  CIRA  model.  Section  III  addresses  the  problem  of  temperature- 
profile  representation.  Section  IV  addresses  the  computational  changes 
needed  to  increase  the  lower  boundary  to  125  km.  Section  V  examines  the 
accuracy  and  efficiency  of  CAD  in  comparison  with  ADEN,  and  Section  VI 


provides  a  summary. 


II.  CIRA  1972  MODEL 


The  CIRA  1972  atmospheric  model  attempts  to  deduce  the  thermal  and 
compositional  structure  of  the  upper  atmosphere  from  the  observation  of 
drag  effects  upon  earth-orbiting  satellites.  This  model  formulates  these 
structures  in  terms  of  a  number  of  variable  parameters,  and  the  values  of 
these  parameters  are  adjusted  to  provide  the  closest  approximation  to  the 
actually-observed  drag  values.  While  it  is  well-recognized  that  this 
solution  is  not  unique  as  far  as  the  thermal  and  compositional  results  are 
concerned  (e.g. ,  several  different  composition  combinations  can  be  derived 
to  produce  good  matches  with  the  drag  data) ,  the  total  mass  density,  being 
directly  related  to  the  measured  drag,  is  represented  rather  well. 

The  CIRA  assumes  a  fixed  temperature  and  fixed  densities  for  the  at¬ 
mospheric  gas  species  at  a  lower  boundary  of  90  km.  Up  to  100  km  a  condi¬ 
tion  of  mixing  is  assumed  to  exist.  Above  100  km  each  constituent  is  as¬ 
sumed  to  be  in  separative  diffusive  equilibrium  such  that  its  number  density 
n  at  any  altitude  z  may  be  obtained  by  integrating  the  diffusion  equation 

din  n  =  -(l+a)dln  T  -  (mg/kT)dz  (1) 

where  a  is  the  thermal-diffusion  coefficient,  m  is  the  species  mass,  g  is 
the  acceleration  due  to  gravity,  k  is  Boltzmann's  constant,  and  T  is  tempe¬ 
rature.  The  integration  of  Equation  1  requires  evaluation  of  the  integral 

I  -fj*  dz  (2) 


where  the  formula 


g  =  9.80665/(1  +  z/R)2 


(3) 


with  R  =  6356.766  km  being  the  radius  of  the  earth,  is  used  to  represent  the 
variation  of  gravitational  acceleration  with  altitude.  The  temperature  pro¬ 
file  formulas  used  by  the  CIRA,  however,  do  not  yield  an  integral  which  can 
be  evaluated  analytically.  This  is  the  primary  task  of  the  present  study: 
to  find  an  alternate  representation  for  the  temperature  profile  which  remains 
a  good  approximation  to  the  CIRA  while  yielding  an  analytical ly-integrable 
form  in  the  diffusion  equation. 

III.  TEMPERATURE  PROFILE 

In  the  altitude  region  above  125  km  the  CIRA  uses  the  formula 


T 

z 

{.0852718 


T125  +  ir(T»  "  T125} 


[(T125-183)/(T,-T125)] 


<*>  125- 


tan 

(z-125)  [l+.  0000045  (z-125)2*5]} 


(4) 


for  the  temperature  profile.  ,  the  temperature  at  altitude  z,  is  asympto¬ 
tic  to  an  "exospheric"  temperature  T^  at  high  altitudes.  T125  9*ven  by 
this  model  in  terms  of  T 

00 

T. __  =  371.6678  +  0.0518806  T 
125  00 

-  294.3505  exp(-0. 00216222  T  )  (5) 
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The  form  of  Tz  in  Equation  4  does  not  yield  an  analytically-integrable  ex¬ 
pression  in  Equation  2.  It  has  been  shown  (Ref.  3,4)  that  a  form  which  does 
yield  an  analytically-integrable  expression  is 

Tz  =  Too  -  (T*.  -  T125)  exp  (-sz  )  -  (6) 


where  s  is  a  constant  "shape  parameter"  and 
z'  =  (z-125)  (R+125)/ (R+z) 


(7) 


is  the  geopotential  height  above  z  =  125  km. 

As  we  are  interested  here  in  altitudes  small  with  respect  to  the  earth's 
radius,  this  equation  closely  approximates  a  simple  exponential-like  asympto¬ 
tic  approach  to  Tot  at  high  altitudes.  The  more  complex  exponential  argument 
is  required  for  exact  integrability  when  the  altitude  variation  of  g  (Equa¬ 
tion  3)  is  taken  into  account.  If  Equation  6  is  adopted  for  the  temperature 
profile,  Equation  2  may  be  analytically  integrated  to  yield  as  a  solution  to 
Equation  1  for  gas  species  i 
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We  thus  find  that  it  will  be  possible  to  replace  the  numerical  integra¬ 
tion  of  the  diffusion  equation  as  employed  by  the  CIRA  code  ADEN  with  a 
simple  closed,  analytic  expression  (Equation  8)  if  Equation  6  is  capable  of 
fitting  the  CIRA  temperature  profiles  sufficiently  accurately.  That  this  is 
possible  has  been  verified  by  fitting  these  CIRA  temperature  profiles  (using 
values  at  25  km  intervals  from  125  to  500  km  altitude)  by  Equation  6  for 
T^,  =  500  to  1900  K  in  100-K  steps,  requiring  T^  and  T12,.  to  be  the  same  for 
the  two  profile  forms  and  finding  the  value  of  s  for  each  value  of  T^  to 
give  the  closest  agreement  between  the  two  profiles.  The  results  of  this 
exercise  are  shown  in  Figure  1.  Here  are  shown  the  profile  results,  with 
the  continuous  curves  representing  the  CIRA  model  and  the  discrete  points 
representing  the  fits.  From  measurements  of  exospheric  temperature  at 
Millstone  Hill  by  the  incoherent  scatter  radar  method  (Ref.  5)  it  has  been 
shown  that  the  random  (unmodelable)  variability  of  T^  is  on  the  order  of 
50  K.  Thus  the  agreements  shown  in  Figure  1  are  considered  acceptable  in 
light  of  the  uncertainty  in  our  knowledge  of  the  instantaneous  state  of  the 
upper  atmosphere. 

Figure  2  plots  s  vs  T^  and  fits  this  trend  by  the  equation 

s  =  ai  +  a2/(Too+a3)  km  1  (9) 

with  the  results  shown  in  Table  1  for  a^,  a2,  and  a3.  The  continuous  curve 
in  this  figure  represents  the  s  values  determined  to  fit  the  CIRA  profiles 
while  the  data  points  represent  the  fit  of  Equation  9  to  these  values. 
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ALTITUDE  (km) 


200  500  800  1100  1400  1700  2000 


TEMPERATURE  (K) 

Fig.  1.  CIRA  temperature  profiles  (solid  lines)  and 
the  fits  afforded  by  Equation  6  (symbols). 


TABLE  I 


FITTED  S  AND  BOUNDARY-CONCENTRATION  PARAMETERS  (EQUATION  10) 


ai 

a2 

a3 

s 

.0057 

17.6 

244 

N2 

17.419 

-110 

133 

°2 

16.583 

-155 

152 

0 

16.944 

140 

3163 

H 

A 

13.397 

97 

335 
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IV.  BOUNDARY  CONDITIONS 

The  CXRA  incorporates  fixed,  non-varying  concentrations  for  the  gas 
species  at  a  lower  boundary  of  90  km.  Because  the  temperature  profile 
in  the  90-125  km  region  varies  according  to  the  solar /geophysical  conditions, 
however,  we  cannot  assume  fixed,  non-varying  concentrations  for  the  gas 
species  at  the  125-km  boundary  adopted  for  the  present  purposes.  Fortunately, 
the  entire  temperature  profile  in  the  CIRA  is  defined  once  T^  has  been 
chosen  such  that  it  is  possible  to  determine  the  concentrations  of  the 
gas  species  at  125  km  altitude  solely  as  a  function  of  T^.  Once  these 
"boundary  variations"  are  obtained,  the  region  below  125  km  may  be  abandoned . 

The  125-km  boundary  variations  for  the  four  CIRA  gas  species  under 
consideration  are  shown  in  Figure  3  as  a  function  of  T  (continuous  curves) . 
For  each  species  these  boundary  "data"  have  been  fitted  by  the  same 
formula  (Equation  9)  used  to  model  the  s  variation.  The  data  points  in 
Figure  3  represent  the  fits  to  the  model  variations.  Table  1  also  lists 
the  results  for  a^,  a2»  and  a3  for  each  gas  species. 
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RELATIVE  CONCENTRATION 


{ 


EXOSPHERIC  TEMPERATURE  OO 


Fig.  3.  Variations  of  the  logarithms  (base  10)  of  the  125-km 
concentrations  versus  exospheric  temperature  for  the  CIRA  model 
(solid  lines)  and  the  fits  afforded  by  Equation  9  (symbols). 
Each  variation  is  normalized  relative  to  its  value  at  500  K. 
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V.  ACCURACY  AND  EFFICIENCY 

The  accuracy  of  CAD  relative  to  the  official  CIRA  code  ADEN  has  been 

checked  by  comparing  the  outputs  of  the  two  routines  for  altitudes  from 

125  to  500  km  (25-km  steps)  and  solar-flux  values  from  70  to  280  units 
-22  -2  -1 

(1  unit  =10  Wm  Hz  )  (10-unit  steps)  with  all  other  inputs  set  to  zero. 
For  no  combination  of  these  height-flux  variable  did  CAD  differ  by  more 
than  7%  from  ADEN.  Below  a  flux  value  of  250  units,  a  4%  error  was  the 
maximum  encountered.  The  average  difference  between  CAD  and  ADEN  was 
about  2%.  This  is  relatively  insignificant  with  respect  to  the  uncertainties 
in  our  knowledge  of  the  state  of  the  atmosphere. 

The  execution  of  one  call  to  CAD  required  4  to  5  ms  on  a  Harris  125 
machine  while  ADEN  required  137  to  140  ms.  Thus  CAD  is  some  30  times 
faster  than  the  CIRA  routine.  ADEN  requires  essentially  the  same  execution 
time  for  any  altitude  from  100  to  500  km.  This  is  because  the  hydrogen 
concentration,  though  of  little  significance  to  total  density  below  500  km 
(it  is  not  even  included  in  the  CIRA  tables  below  this  altitude) ,  is  included 
in  the  ADEN  calculation  to  avoid  any  discontinuity  at  this  altitude.  Be¬ 
cause  the  CIRA  adopts  a  boundary  variation  for  the  hydrogen  concentration 
at  the  500-km  level,  the  inclusion  of  hydrogen  at  lower  altitudes  requires 
the  integration  of  the  diffusion  equation  downward  from  500  km  to  the 
altitude  of  interest.  Thus  the  numerical  integration  of  ADEN  must  cover 
the  entire  100-to-500-km  range  no  matter  which  altitude  in  this  region  is 
under  consideration. 
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VI.  SUMMARY 


Procedure  CAD  is  intended  to  be  an  efficient  machine  implementation 
of  the  CXRA  1972  atmospheric  drag  model,  designed  specifically  for  use 
in  the  approximation  of  atmospheric  drag  effects  upon  orbiting  satellites. 
The  code  ADEN  provided  by  the  CIRA  publication  for  this  calculation  is 
altered  to  (a)  approximate  the  numerical  integration  of  the  diffusion 
equation  by  an  analytic  expression,  (b)  limit  the  altitude  range  of 
calculation  to  the  satellite-orbit  region  in  which  significant  drag  effects 
are  experienced,  and  (c)  eliminate  minor  constituents  of  insignificant 
density  contributions.  These  modifications  render  CAD  some  30  times 
faster  them  ADEN  while  accuracy  is  maintained  to  within  a  few  percent. 

The  appendix  presents  a  FORTRAN  implementation  of  procedure  CAD, 
retaining  the  form  of  ADEN  as  much  as  possible.  The  input  and  output 
arguments  for  the  two  routines  are  identical. 
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SUBROUTINE  CAD  (CALCULATE  ATMOSPHERIC  DESITY)  IS  A  MODIFICATION  OF 
SUBROUTINE  ADEN t  UHICH  IS  THE  FORTRAN  PACKAGE  INCLUDED  IN  THE  CIRA  (COSPAR 
INTERNATIONAL  REFERENCE  ATMOSPHERE)  1972  PUBLICATION  (ACADEMIE  VERLAG » 
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